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ABSTRACT 



This paper presents an accuracy analysis of a suggested approximate 
confidence interval for system maintainability parameters. 

Technically, the simulation demonstrates feasible ranges of parameter 
applicability for a fit of linear combinations of generated gamma variates 
to the gamma distribution, using the method of moments. 

The simulation has application to the classical confidence interval 
for mean time to repair of a series system^ under the assumptions of 
gamma distributed repair times, and method of moments estimators,. 

The paper provides no validated conclusions although it does display 
parameters and ranges of apparent extremely high model validity. 
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I . INTRODUCTION 



The purpose of this paper is to present a computer simuTatioir. in 
order to investigate the adequacy of a suggested approximate confxdence 
interval for system maintainability parameters. In a technical sense, 
it will demonstrate feasible ranges of parameter applicability fbn a fit 
of linear combinations of generated gamma variates to the gamma, distribu- 
tion, using the method of moments. The application of the simulation is 
to the classical confidence interval for mean time to repair' a system, 
under the assumptions of series components, gamma distributed repair 
times, and method of moments estimators. The parameters of the distribu- 
tions that Mere investigated were selected so that the maintainability 
curves would be representative of those curves presently of interest in 
Industry, 
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II. SUMMARY 



A procedure for a fit of linear combinations of generated gamma 
variates to the gamma distribution, using the method of moments is 
presented in Chapter III. A simulation of this procedure is proposed 
in Chapter IV and the results are tabulated in Tables I and II of' 
Chapter V. 
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III. THE STATISTICAI. MODEL 



A. SERIES ASSUMPTION 

Suppose we are given a system of components which are not necessarily 
series, where the time to failure of the ith component of the system has the 
exponential (X^) distribution. 

Then , 

1) for the case where there are identical components in. series of 
type i within the system, set X^^ = 

2) for the case where there are h^ identical components in parallel 
of type j within the system, set X^ ' = 

Therefore, for some parallel-series combination there is obtainable 
a series system of uniquely different components whose failure ratios are 
assumed to be exponential. 

Thus, the system being described is a series system with k different 

components and failure rates X , X , •••, X.**, X/, \ where 

1 z X j k. 

all the X. are in the same units. 

1 



B. FAILURE RATE ASSUMPTION 

Let denote the time to failure of component i, where has the 
exponential (X^) distribution. Thus, 

f (x. ; X.) = X.e ^i^i (3-1) 

X. 1 X X 

X 

We further assume that the system fails when exactly one component 
fails or the component which caused the system failure takes the longest 
time to repair. 
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C. GAMMA ASSUMPTION 



Let denote the time required to repair the ith component and 
suppose that has the gamma distribution. Thus, "" T a.. ,, h.).. 

Therefore, 



b. b.-l -a.t. 
a. 1 t . 1 e 1 r 

^T. " r(bj 

1 1 



(a-2) 



is the density equation with 



E[T.] = = y. 

1 a. 1 
1 



(a-3) 



which is the mean time to repair (MTTR) for the ith component. 

D, SYSTEM MTTR OBTAINED 

Let 0 denote the mean time to repair for the system. Tlierefore,, 



9 = Z P [component i fails first] u. 
i ^ 



( 3 - 4 ) 



and, if the assumption of a series system is valid, then 



X = Z X. 



i=l 



and 



k X. b. 

0 = MTTR = E (-r^) (— ) 
. 1 X a . 

1=1 1 



(3-5) 



E. CONFIDENCE LIMIT OBTAINED 

An estimated upper confidence limit for system maintainability, denoted 

by 9^, has been proposed. 9^ was derived by the following procedure: 

For any component i , T_ _ , T. _ , T. ^ , ••• T. is a random sample on T., 

ii ±Z iJ in. 1 

1 

the time to repair component i* Failure data is also available on component 



i, so that an estimate X. of X. is possible. 
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Thus, 9 can be estimated by 9 where 



k X. _ 
9 = Z T. 
1=1 * 



(3-6) 



and 

T. = — Z T. . 

1 n . . 11 

1 J 

Moreover, because of the large amount of rndustrial testing,, it:, can 
be assumed that = X^, and for purposes of derivation we shall, treat 
the X^ as though they are constants • Thus, 



k X. _ 
9 = Z — T 
i=l ^ 



(3-7) 



k X. b. 

E [9] = Z ^ (-5^) 
i=l ^ 



(3-8) 



and 



k X. 2 b. 

var (9) = Z (^) (^) . (3-9) 

i=l a. 

1 

We shall fit a two parameter gamma to the density of 0 and obtain 
the upper confidence limit from this fitted distribution. Thus we are 
assuming 0 '' r(a,b) and use will be made of the method of moments to make 
the fit. Through use of formulas (3-8) and (3-9) 



k X. b. 

E [9] = S (t^) = 7 (3-10) 

, - A a . a 
1=1 1 



and 



var (9) 



k 

Z 

i=l 



X. 2 b., 

C^) (^) 

a^ 

1 



h 

Z 

a 



( 3 - 11 ) 
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The parameters a and b can be solved for, as follows: 



a 



k 

E 

i=l 

k 

E 

i=l 



X. b. 

<r> 

1 



X. b. 

a . 

1 



(a-iz) 



and 



b = 



k X . b . 

^ ir ~ 

1=1 ^ 

k X. 2 b. 
1=1 ^ a/ 



For the available data, this then becomes: 



Cl^'3) 



where , 



k 

E 



k 

E 

1=1 




2 




S. 

1 



n. 



1 

E 



-1 “ 

"l ^ j=l ^ 




(3-14) 



2 

If b is an integer, 20a has the X^2b) therefore,, define 

(2b) as the integer closest to 2b and the following approximate confidence 
interval can be formed: 
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1 - a = P 



C3-15) 



’ X(l-a)(2b)’ 

P r20^ ^ ^(l-g)(2b) 

^(2b) (2b) J 



= P 



b ^ 9 (2b) 

2 

^(1-a) (2b) 



= P 



fb 9 (2b) 

2 

^(1-a) (2b) 



C3-16) 



where the approximate equality compensates for the use of estimators. 
Finally, the choice of a = .20 results in the following: 




9 (2b) 

2 

^(.80)(2b) 



C3-I7) 



Thus, the desired suggested approximate 80% confidence limit for system 

" " 2 

maintainability is 0(2b)/x/ o\/nZ\- 

C Czb; 
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IV, THE SIMULATION PROCEDURE 



As explained in Chapter III, the system to be simulated cons±sts of 
k components in logical series with system maintainability (&) expressed 
as : 



k X. b. 

» - <r> 

1=1 1 



(^- 1 ) 



where is the true maintainability of the ith component of the 

system. 

Denote an upper confidence limit for 0 by 0 . If 0 is in. fact the 



u 



exact 100 (l-a)% upper confidence limit for 0, then 



P 




0 (2b) , ^ 

2 J ^ 

^(l-a)(2b) 



a 



holds to a reasonably close approximation* 

In fact, b/a should then be the ath percentile point of the simulated 

distribution of 9 * 
u 




b 

a 



The choice of a = .20 defined G as the 20th percentile point of 

u(2) 

the distribution of 9 . 

u 
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In order to investigate this distribution, a computer was utilized 
to generate the required gamma variates* Then, 500 values of 0- were 

XX 

computer and ordered such that 



®ul " ®u2 " ®u3 < 



< 9 



C'4^Z) 



u500 

Since it was desired to display the 80% upper confidence limit for 

0 (0 /o \)5 which implies 80% of the ordered values be greater than 9 

u(.2; ^ ° n(.2) 

the 20th percentile point of the 9^ distribution was found. This ZQth. 
percentile point is the 100th ordered value in formula (4-2) above and, 
if the procedure is correct, should approximately equal 

a 

The primary measure of accuracy for the simulation will be the value 
of 
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u lOO a ' 
b/a 



(4-3) 



which is an expression relating the estimated value of system maintain- 
ability to the gamma value for MTTR (b/a)* Thus, the accuracy 

of the simulation is presented in the notion of relative error. 

Also, the statistic 



[#9 . > b/a] 
500 



(4-4) 



will be computed in order to display the relative error between the true 
value of system maintainability (b/a) and the number of generated esti- 
mates of this value (9 .). 

wj 

The analysis for this proposed method was conducted using different 
combinations of values of the gamma input parameters, varying values,, 
one of two values for k, and n^ values of 10, 20, 50, and 100, 
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A, THE SIMULATION PROGRAM EXPLAINED 



Available in Appendix E is a flow chart for the computer simulation.. 
An explanation of the blocks on that diagram follow: 

1) Dimension the computer matrix as required to include the poHSibXe 
values for the following input parameters: 
k 

^ 1 ’ ^ 2 ’ \ 

(a^ y b , ^^2 ’ ^2^^ ***5 ^^k ’ ^k^ 



n^ , n 



1 > - 2 ^ 



• • • , n . 



2) Generate the following n^ random repair times according to the 
gamma distribution with parameters (a^, b^. ). Utilize the IBM Scientifie 
Subroutine Package GAMMA to get: 



11’ 


T^2’ 


T 

13’ 


... T' 

’ In 




T 


T 


. . . T 


21’ 


22’ 


23’ 


’ 2n 



T T T 

k2^ ^k3’ 





3) For each row above, compute the mean value (T^) and the variance 
) to get the pairs : 





2 

k 



4) Utilize formula (3-14) to compute the method of moments estimate 

A 

of the unknown parameter b, denoted by b. 

2 

5) Compute the integer closest to 2b for use with the X^2b) 
distribution. 
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6) Utilize formula (3-6) to compute the estimate of system maintain- 
ability 9, denoted by 0. 

2 

7) Utilize formula (3-16) and the x formal a (see Appendix C) to 

compute the approximate upper confidence Limit on system maintainability. 

Denote this value as 0 . 

u 

8) Repeat step 2 through step 7 a total of 500 times. 

9) Utilize the IBM SSP SHSORT to order the 500 values of" 0 from 

u: 

smallest to largest. 

10) Pick out the 100th value of the above ordering print 

out this value as the 20th percentile point of the distribution.. 

11) Compute the value of b/a and then utilize formula (4-3) to compute 
the primary measure of accuracy for the simulation. Print out this value. 

12) Utilize formula (4-4) to compute another measure of simulation 
accuracy. Print out this value. 
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V. RESULTS 



For the procedure as presented in Chapter III, the value of' system 
maintainability will follow the gamma distribution. The parameters of 
the distribution were chosen so that the maintainability of' each compon- 
ent was at a preselected value position on the curve or would be repre- 
sentative of those curves currently obtainable in industrial! applications 
(See Appendix A.) The following cases were simulated : 



CASE # 


INPUT PARAMETERS- 


I 


k = 15 


II 


= 10, b. =1 
i 1 

= .005 

= 10, 20, 50, 100 
k = 15 


III 


a. = 10, b. = 2 

1 X 

X. = .005 
1 

n, = 10, 20, 50, 100 

X 

k = 15 


IV 


a, = 5, b. = 1 
X, = .005 

X 

n. = 10, 20, 50, 100 

X 

k = 15 

a. = 30, b. = 1, i = 1, 2, '*•, 10 
a^ = 10, b^ = 1, i = 11, 12, ••*, k 
X. = .005 

X 

n. = 10, 20, 50, 100 

X 
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CASE # 



INPUT PARAMETERS 



V 


k » 15 


VI 


= 30, = 1, i = 1, •••,10 

a,=10,b.=l,i=ll,^^^,k 
= .01 

n, = 10, 20, 50, 100 
1 

k = 30 


VII 


a. = 30, b. = 1, i = 1, •••, 10 

a. = 10, b. = 1, i = 11, ••., k 

X. = .005 
1 

n. =10, 20, 50, 100 
1 

k = 30 


VIII 


a^ = 30 , b^ = 1 , i = 1 , • • • , 10 

a. = 10, b. =1, i = 11, •• • , k 
1 1 

X. = .01 
1 

n. =10, 20, 50, 100 
1 

k = 15 

a. = 5, b. = 1 
1 1 

X. = .05 
1 

n. =10, 20, 50, 100 
1 


Each case 


is specified by the relevent input parameters and the vary- 


ing number of 


random repair times generated (n^) . The cases differ 



most significantly in their values for the parameters of the Gamma 

distribution. However, the value of X. is varied as is the value of 

1 

k, the number of components in series. Of course, all of the cases were 
run over an identical range for 
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The results are compiled in the following tables. The column 
headings are self-explanatory and their use has been discussed earlier 
in this report. 
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RELEVENT STATISTICS OF COMPUTER SIMULATION FOR CASES 
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RELEVENT STATISTICS OF COMPUTER SIMULATION FOR CASES V - VIII 
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APPENDIX A 



THE GAMMA DISTRIBUTION 

Use of the gamma function to describe mat rrtenance— time functions 
has been described in the literature. 

The gamma distribution is a two-parameter function. 





a.,bi) 



b. 

1 t. 



h.-l 



t r 



r(b.) 



-a. t . 

X X 

e 



The parameters are denoted by the letters a and b,. Of the two,.b is 
considered to be the more critical because it controls the shape odf the 
curve; while a merely determines the scale of the axes. 




Figure 1. Gamma density function; 

a = 1, b = 0, 1^ 3, 5 

For the distribution as presented, the parameters a and b am 
restricted by the inequalities a > 0, b^Q. Also,, both a and b: will 
remain as whole numbers. 
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The b parameter is seen to vary as a function of the repair concept. 
When the repair concept being described invalves relatively simple 
actions (e.g., black box replacement^ etc..) and very few long: downtime 
tasks, the b parameters will be a small nnmbeir.. As the downtime- density 
involves more time-consuming actions^ tbe b p^arameter. wilX take onr a 
larger value. 

Throughout the literature it has been shown that the gamma": distribu- 
tion could be used to describe a family of distribution curves varying 
from the negative exponential to the normal.. Thus, by suitable choice 
of the parameters a and b, the downtime of a wide variety of:: reasonably 
complex equipments can be described.. 
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APPENDIX B 



THE EXPONENTIAL DISTRIBUTION 

One of the most widely used distributions in fields of Reliability/ 
Maintainability is the one-parameter exponential function*. 

The function is defined by: 



-X.X 

f (X. ,X.)=X.e^^ X>0 

1 X X — 

= 0 elsewhere 

where is the failure rate associated with component i.. 

X. 

-X.t. -X.X. 

F(X.) = P(X < X.) = X. 1 e ^ ^ dt. = 1-e ^ ^ X. > 0 

X — X X J X r* 

° = Q 0 



is the cumulative distribution functiou, and 



E[X.] = 1/X.. 

1 X 

Thus, the expected value equals the reciprocal of the parameter X^. 
var (X.) = 1/X.^ 

X X 




Figure 2. The exponential distribution 
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NOTES : 



1) The exponential distribution is a special case of both the gamma and 
weibull distributions. 

2) It is the distribution which is expected when the mechanisms are sa 
complex that many deteriorations with different failure rates are 
operable . 

3) When parts have an exponential failure distribution, the equipment: 
consisting of these parts also has an exponential distribution function. 
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APPENDIX C 



THE CHI-SQUARE DISTRIBUTION 

A. quite useful case of the gamma distribution occurs when a = T/Z 
and b = n/2, where n is a positive integer. Then, a one— parameter* 
family of distributions is obtainable with density function: 



f (z) 



(n/2 - 1) -z/2 

z e 

2 ”^^ r (n/2) 



z > Q 



A random variable z having the above density is said to have the 

2 

chi-square distribution with n degrees of freedom (denoted by 

In Figure 3 below, the density function for n = 1, 2, and n > 2 is 
shown. 



f U) 



f (z) 






z 



Figure 3 

with E[z] = n, var (z) = 2n. 

Our interest in the chi-square distribution is based upon its many 
important applications in statistical inference. 
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The chi square distribution is tabled for degrees of freedom up to 

30. Thus we may find in the table that value, denoted by ^ ^ and 

^ (a)(n) 

2 

satisfying P(Z y^) = a, 0 < a < 1. 



f (z) 




Figure 4. The chi-square distribution 

For the case where the degrees of freedom exceed 30 (n > 30) the 
chi-square distribution can accurately be approximated by the normal 
distribution as indicated in the following theorem: 

2 

Theorem: Suppose that the random variable Y has distribution y • 

n 

Then for sufficiently large n, the random variable yp2X has approximately 
the distribution N(/2n - 1, 1). Therefore, 
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P(Y < t) = P(v^ < /^) 



= V(/2Y - /2n - 1 ± - /2n 

= i(/2t - /2n - 1) 

where '!> values are obtained from the normal tables. 
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APPENDIX D 



THE METHOD OF MOMENTS 

The oldest general method for generating estimates of unknown: 
parameters, given a sample, is known as the method of moments. Dr is 
a comparatively simple method and it generates "reasonable" estimators. 

In general, the method proceeds as follows: 

Given a random variable X, with distribution function F where, this 
distribution is indexed by the unknown parameter X. Assume the_ first 
moment of X (its mean value) is dependent in some simple way upon X, 
like = g(X). Then, given a sample of n values of X, define the first 
sample moment as X. The method of moments now specifies that X be equated 
to g(X) and finally to solve for X. This resulting value for X is the 
method of moments estimator of X. 

For the two dimensional case, such as X being distributed according 
to the gamma distribution with unknown parameters a and B ,. there is 
required a multi-dimensional parameter space which in the gamma 

case has k = 2. Then, this space can be defined such that 
Q = {(a, B); a > 0, B > 0}cs and further, it can be shown that 

2 2 2 2 
li = aB and y = aB + a B 
X X 

2 

where y is the second moment with respect to X. Then,, as in the one 
A 

dimensional case, a random sample is required which can be analyzed to 
give : 

X = aB 

2 —2 2 2 2 
+ X = aB + a B 

A 
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which allows , as shown in Chapter III 



a 





as acceptable estimators of the parameters a and 6. 
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APPENDIX E 



© 



READ IN 
PARAMETERS 



SUBROUTINE 

GAMMA 



GENERATE 

GAMMA 

VARIATES 



COMPUTE 
T and s' 



COMPUTE 
b and 2b 




SUBROUTINE V 

SHSORT / 



ORDER 500 
VALUES OF 

u(i) 



PRINT 20th 
PERCENTILE 
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oo oooo non noon o on 



THIS PROGRAM COMPUTES THE 80 PERCENT UPPER CONFIDENCE LIMIT 

DIMENSION THE MATRIX 

P.EAL^A LAMBDA! 100) ,N( 100) ,LAM 

DIMENSION T(ino,100) ,Te(100),S2(100) , TUI 500) , ALPHA! 100) , BETA ! 10 

$0),KEY!500) 

READ INTO PROGRAM PARAMETERS OF INTEREST 
1235 REA^!5,lno,E^'^=4321 ) K, (LAMBDA! I ) ,I = 1,K) , I ALPHA ! J ) r BETA ( J ) ,-J= 1 , K ) 
REA0!5,105) NI 

PROGRAM HILL LOOP 500 TiMESt GENERATING INDEPENDENT ESTIMATES OF 
THE 80 PERCENT UPPER CONFIDENCE LIMIT 

DO 1234 JJ=1,500 

GENERATE GAMMA TIMES TO REPAIR 

DO 1 1=1, K 
N! I ) = NI 
NN=N( I ) 

DO 1 J = 1 , NN 

CALL GAMMA! BETA! I ) , ALPHA! I ) ,X ) 

T! I ,J) = X 

1 CONTINUE 

COMPUTE MEAN TIME TO REPAIR FOR EACH K GAMMAS GENERATED AMD THE 
VARIANCE FROM THE METHOD OF MOMENTS FIT 

DO 2 I=1,K 
SUM=0.0 
SUMS0=0.0 
NN=N!I ) 

DO 3 J=1,NN 
SUM=SUM + T!I,J) 

3 SUMSO=SUMSQ + T ( I , J ) *T ( I , J ) 

TB( I ) = SUM/N! I ) 

S2(I) = (SUMSQ-N! I )*TB( D^TB! I) )/!N!I )-1.0) 

2 CONTINUE 
LAM=0.0 

DO 4 1=1,K 

LAM= LAM+LAMBDAI I ) 

4 CONTINUE 



SUM=0.0 
SUMSC=0.0 
DO 5 I =1 ,K 

SUM= SUM + (LAMBDA! I) /LAM):i=TB! I ) 

SUMSQ = SUMSQ + ! LAMBD A! I ) / LAM ) =<= ! LAMBDA ( I ) /L AM ) «( S2 ( I ) *S2 ! I ) ) 



31 



on noon oooo on non ooo ooo ooo 



5 CONTINUE 

COMPUTE VALUE FOR ESTIMATE OF PARAMETER BETAr CALLED BHAT 
8HAT=(SUM*SUM) /SUMSQ 

COMPUTE INTEGER CLOSEST TO TWO BHAT 

IH=2.0*BHAT+0.5 
BHAT2 = IH 

COMPUTE ESTIMATE OF MEAN TIME TO REPAIR FOR SVSTEMtCALL IT THETA 

THETA = 0,0 
DO 6 1=1 ,K 

THETA = THETA + ( L AMB0A( I )/LAM I=f=TBf I ) 

6 CONTINUE 

UTILIZE CHISQ FORMULATION TO GET STATISTIC FOR CONFIDENCE LIMITS 
CHISA=(2.0/(9.0^BHAT2) ) 

CH I S0=BHAT2=!' ( 1 . C-CHI SA-C . 341 78*SQRT ( CHI SA ) 

TU( JJ) = ( THETA=:=BHAT2)/CHISQ 



1234 CONTINUE 

NOW THAT THE 500 ESTIMATED VALUES HAVE BEEN GEfJERATED, SORT THEM 
INTO ORDER 

KK=500 

DO 99 1=] ,KK 
99 KEY( n = I 

CALL SHSORT(TU,KEYtKK) 

PICK OUT AMD PRINT THE 2CTH PERCENTILE ACCORDING TO THE FOLLOW- 
ING FORMAT 

WRITE(6,222) M I , Kt LA'^BDA ( 1 ) 

WRITE(6,200) TU(IOO) 



222 FORMAT (• 1» , 'Ml = »I5, ’ K=*tI5, ' LAMSDA= • , F 10 . 5, // / ) 

100 F0RMAT(I5/ 15F5.3 / 15F5.3 / ( 1 0 ( 4. 0 , F4 . 0 ) ) ) 

105 FORMAT ( 15 ) 

200 FORMAT (IX, • THE 20TH PERCENTILE IS ',F15.5) 
WRITE(6,400MTU( I ) , I = I ,500) 
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400 FORMATC • , 10 ( F 8 .4 , 2X ) ) 

C 

C COMPUTE THE VALUE OF THE THEORETICAL MEAN, CALL IT BOA 

C 

DEN0H=0.0 
TOP=0.0 
00 50 1=1, K 

TOP = TOP +(LAMBDA(I)/LAM)=!=(BETA(I )/ALPHA(I) ) 

OENOM = nENOM+(L A^BDA ( I ) /L AM ) -( 6ETA( I ) / ALPHA! I ) *-2) 

50 CONTINUE 

A = TOP /D£ MOM 
B = T0P*-'r2/DEN0M 
BOA = B/A 

C COMPUTE THE ACTUAL PERCENTAGE 

C=0.0 

DO 60 1=1,500 
IF (TU( I ) .GE.BOA) C=C+1.0 
60 CONTINUE 

STAT = C/500.0 

C COMPUTE THE ABSOLUTE VALUE OF THE RELATIVE DIFFERENCE 

E= ABS(TU( L00)-B0A)/B0A 

C PRINT our the ABOVE STATISTICS , ACCORDING TO THE FOLLOWING FORMAT 

WRITE(6,220) BOA, STAT, 3, A 

220 FORMATdX,’ B OVER A IS ',F15,5,5X,' STAT IS ',F15.5r5X,' B IS 
$F15.5,5X,» A IS *,F15.5) 

WRITE(6,230) E 

230 FORMATdX,' ABS OF VALUE IS ',F15.5) 

IN = 20 

RANGEd ) = 0.1822 + 0.00159 

DO 91 1=2,20 

91 RANGEd) = RANGEd-1) + 0.00159 
K=1 

DO 98 1=1,20 
FREQ( I) = 0.0 
DO 97 J=K,500 

IF(RANGE( T ) .LT .TU( J ) ) GO TO 96 
FREQ( I )=FREQ( I )+1.0 
GO TO 97 

96 K=J 

GO TO 98 

97 CONTINUE 

98 CONTINUE 
GO TO 1235 

4321 STOP 
END 
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SUBROUTIME GAMMA(BtA,X) 
K=B 

TR=1.0 
DO 5 1=1, K 
R=RN(0) 

5 TR=TR*R 

X=-ALOG(TR)/A 

RETURN 

END 



:C= ^ 5k sjc ^ ^ ^ ^ * 



PURPOSE 

RN COMPUTES UNIFORMLY DISTRIBUTED RANDOM VARIATES OVER THE RANGE 
(0,1). THE PERIOD FOR THIS GENERATOR IS 2=!=^32 WHICH IS THE LENGTH 
OF THE FULL SY5TEM/360 WORD GIVING TWICE THE PERIOD OF RANDU. 

USAGE 

IN X = RN(0) WHERE X WILL BE ASSIGNED THE NEXT RANDOM DEVIATE. 

DESCRIPTION OP PARAMETERS 
NO FORMAL PARAMETERS 

REMARKS 

THIS ASSE’‘*B!.Y LANGUAGE ROUTINE GIVES TWICE THE PERIOD OF EXISTING 
GENERATORS AND IS EXTREMELY FAST. 11.88 MICROSECONDS WITH 4.23 
MICROSECONDS FOR THP CALLING SEOUENCE. IT IS TAKEN DIRECTLY FROM 
, COMM. ACM 12,12 (DEC. 1969) 695. 



METHOD 

LEHMER’S MULTIPLICATIVE CONGRUENTIAL SCHEME IS USED. 

U(K+1) = L*U(K) (MOD P), X(K+1) = U(K+1)/P WHERE P = 2**32 AND 

L = 32781. 



RN 



MULT 

ZERO 

FLOT 

SEED 
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CSECT 

USING 

L 

M 

ST 

LD 

AD 

BR 

DC 

DC 

DC 

DC 

END 



*,15 

1 , SEED 

0,MULT 

I , SEED 

.0,FL0T 

0,ZER0 

14 

F'32731 * 

D’O* 

X'46000000 

F* 1* 

RN 



SET BASE ADDRESS 
LOAD U(K) 
SEED=SFED*MULT 
STORE U(K+1) 
FLOAT AND 
NORMALIZE 
RETURN 
2**32 



RN 
RN 
RN 
RN 
R N 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RM 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RN 
RM 
RN 
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